function obj=exp_u_delta(y, rho, kappa, zeta, logit_coef, X, delta)


X_new=[y'; X]; 

temp=prod(y.^zeta,2);

prob=exp(sum(logit_coef.*X_new,1))/(1+exp(sum(logit_coef.*X_new,1)))*delta;

obj=prob*u(rho,temp-kappa)+(1-prob)*u(rho, temp);


end